# Replication code for Daniel Lobo and Ryan Brutger's Appendix (Last updated October 20, 2024)
# "Fairness According to Whom"
# Code for the manuscript is in the companion R script


# Code was run on R version 4.0.1 (2020-06-06) -- "See Things Now"
# Using MAcBook Pro using macOS Catalina, with 2.9GHz Dual-Core Intel Core i5

# Race/Ethnicity Codings
# 1= White
# 2= African American
# 3 = Hispanic
# 4 = Asian
# 5 = Native American
# 6 = Pacific Islander
# 7 = Other

################# 1. Load libraries, load and clean data ##############
# Please note, you will need to run the "interaction_plots.R" (by Anton Strezhnev) prior to generating the marginal effects plot

#First, set your working directory using the setwd() command:
setwd()

#Load data
data <- get(load("Fairness_to_Whom_Replication.RData"))

#Load required R packages (download first if necessary)
library(foreign)
library(Hmisc)
library(ggplot2)
library(psych)
library(stargazer)
library(mediation)

# Create indicator for white respondents
data$white <- ifelse(data$race==1, 1, 0)
# Sreate indicator for men
data$men <- ifelse(data$gender==1, 1, 0)

# Subset to only white and Black respondents
data <- subset(data[data$race<3, ])

# Demographics for the appendix for the "meaning of fairness" study
# Table 1 of appendix

# Create income brackets
data$Income_0_50 <- ifelse(data$income==1 | data$income==2, 1, 0)
data$Income_50_100 <- ifelse(data$income==3 | data$income==4, 1, 0)
data$Income_100_150 <- ifelse(data$income==5 | data$income==6, 1, 0)
data$Income_150up <- ifelse(data$income==7 | data$income==8 |data$income== 9, 1, 0)

# Create party variables
data$Democrat <- ifelse(data$party==1, 1, 0)
data$Republican <- ifelse(data$party==2, 1, 0)
data$Independent <- ifelse(data$party==3, 1, 0)

# Create age brackets
data$Age18_29 <- ifelse(data$age<30 & data$age>17, 1, 0 )
data$Age30_44 <- ifelse(data$age<45 & data$age>29, 1, 0 )
data$Age45_59 <- ifelse(data$age<60 & data$age>44, 1, 0 )
data$Age60up <- ifelse(data$age>59, 1, 0 )

# Age variable
table(data$age)
data$age0 <- ifelse(data$age>= 18 & data$age <= 29,1,0)
data$age1 <- ifelse(data$age>= 30 & data$age <= 44,1,0)
data$age2 <- ifelse(data$age >= 45 & data$age <= 64,1,0)
data$age3 <- ifelse(data$age >=65,1,0)
data$agerange <- ifelse(data$age0==1, 1,ifelse(data$age1==1, 2,ifelse(data$age2==1, 3, ifelse(data$age3==1, 4, NA))))

# Race variables
# Create indicator for white respondents
data$White <- ifelse(data$race==1, 1, 0)
data$Black <- ifelse(data$race==2, 1, 0)

#  subset to those who answered the "what does fairness mean" question
dataF <- subset(data, is.na(data$fair4)==FALSE )

# Produce .tex code for Table 1 of appendix, Meaning of Fairness Study Demographics
stargazer(dataF[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women", "White", "Black", "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)

# Analyze the Trade Concessions Study
# "V1" is used to indicate the data is from the Trade Concessions study

# Create data subset to those who answered the support dependent variable
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE ) 
# Create data subset those who answered fairness and support question
dataV1_fair <- subset(dataV1[is.na(dataV1$V1_fair)==FALSE, ])

# V1 - subset to those who answered  fairness dependent variable
dataV1 <- subset(data, is.na(data$V1_fair)==FALSE )

# For Table 2 of the appendix, the N is manually added based off the number of respondents in each relevant group
# Produce .tex code for Column 1 of Table 2 of appendix
stargazer(dataV1[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women", "White", "Black", "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)

# V1 Limit to Favorable Condition
dataFav <- subset(dataV1, dataV1$V1_fav==1 )
# Produce output for column 2, Favorable Cond.
stargazer(dataFav[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women", "White", "Black", "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)
# V1 Limit to Unfavorable Condition
dataUnFav <- subset(dataV1, dataV1$V1_unfav==1 )
# Produce output for column 2, Favorable Cond.
stargazer(dataUnFav[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women", "White", "Black", "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)
# V1 Limit to Equal Condition
dataEqual <- subset(dataV1, dataV1$V1_equal==1 )
# Produce output for column 4, Equal Cond.
stargazer(dataEqual[c("Age18_29", "Age30_44", "Age45_59", "Age60up", "women", "White", "Black", "Income_0_50", "Income_50_100", "Income_100_150", "Income_150up", "Democrat", "Republican", "Independent")], type = "latex", title="Study 1 Demographics", 
          summary.stat = c("mean"), digits = 2)

# Number of Black respondents in each sub-treatment for appendix Table 2.1
# Look at how many white and Black respondents are in each sub-treatment condition

# Subset data by respondent's race
dataV1W_fair <- subset(dataV1_fair, dataV1_fair$race==1 ) 
dataV1B_fair <- subset(dataV1_fair, dataV1_fair$race==2 ) 

# Favorable Treatments
sum(dataV1B_fair$T9030[dataV1B_fair$T9030==1]) # 21
sum(dataV1B_fair$T9060[dataV1B_fair$T9060==1]) # 25
sum(dataV1B_fair$T6030[dataV1B_fair$T6030==1]) # 17 

# Unfavorable Treatments
sum(dataV1B_fair$T3090[dataV1B_fair$T3090==1]) # 19, 9% of Black respondents
sum(dataV1B_fair$T6090[dataV1B_fair$T6090==1]) # 28, 13% of Black respondents
sum(dataV1B_fair$T3060[dataV1B_fair$T3060==1]) # 18 

# Equal Treatments
sum(dataV1B_fair$T3030[dataV1B_fair$T3030==1]) # 31 
sum(dataV1B_fair$T6060[dataV1B_fair$T6060==1]) # 25
sum(dataV1B_fair$T9090[dataV1B_fair$T9090==1]) # 27


# For analysis in section 2.2 of appendix
# looking at fairness evaluations across sub-treatments 

# Set critical value used for confidence intervals
c.value <-qnorm(.95)

# Favorable conditions
mean(dataV1B_fair$V1_fair[dataV1B_fair$T9060==1]) # 0.12, this falls between the means for the other two favorable treatments
mean(dataV1B_fair$V1_fair[dataV1B_fair$T6030==1]) # 0.353
mean(dataV1B_fair$V1_fair[dataV1B_fair$T9030==1]) # -0.238

# Unfavorable conditions
mean(dataV1B_fair$V1_fair[dataV1B_fair$T6090==1]) # -0.036
mean(dataV1B_fair$V1_fair[dataV1B_fair$T3090==1]) # 0.00
mean(dataV1B_fair$V1_fair[dataV1B_fair$T3060==1]) #0.222

# Analysis for section 4 of appendix, Figure 1 of appendix
# Create Figure for men and women's fairness scores by treatment
dataV1men_fair <- subset(dataV1_fair, dataV1_fair$gender==1 ) 
dataV1women_fair <- subset(dataV1_fair, dataV1_fair$gender==2 ) 

# Generate estimates for men
V1_equal <- mean(dataV1men_fair$V1_fair[dataV1men_fair$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1men_fair$V1_fair[dataV1men_fair$V1_equal==1])/length(dataV1men_fair$V1_fair[dataV1men_fair$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1men_fair$V1_fair[dataV1men_fair$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1men_fair$V1_fair[dataV1men_fair$V1_unfav==1])/length(dataV1men_fair$V1_fair[dataV1men_fair$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1men_fair$V1_fair[dataV1men_fair$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1men_fair$V1_fair[dataV1men_fair$V1_fav==1])/length(dataV1men_fair$V1_fair[dataV1men_fair$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate fairness estimates for plot
V1.means.Menf <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Menf <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Menf <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Menf <- c("Equal", "Unfavorable", "Favorable")

# Generate estimates for women
V1_equal <- mean(dataV1women_fair$V1_fair[dataV1women_fair$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1women_fair$V1_fair[dataV1women_fair$V1_equal==1])/length(dataV1women_fair$V1_fair[dataV1women_fair$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1women_fair$V1_fair[dataV1women_fair$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1women_fair$V1_fair[dataV1women_fair$V1_unfav==1])/length(dataV1women_fair$V1_fair[dataV1women_fair$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1women_fair$V1_fair[dataV1women_fair$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1women_fair$V1_fair[dataV1women_fair$V1_fav==1])/length(dataV1women_fair$V1_fair[dataV1women_fair$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate fairness estimates for plot
V1.means.Womenf <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Womenf <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Womenf <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Womenf <- c("Equal", "Unfavorable", "Favorable")

# Figure: Men/Women Fairness scores by Equal, Unfav, Fav treatments
par(mar=c(3,6,1,2))
plot(1:3, V1.means.Menf , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .9), xlim = c(0.75, 3.25), pch = 22, type = "p", xaxt="n", cex=1)
mtext("Average Fairness ", side=2, at=.25, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.Menf)){
  lines(c(i,i), c(V1.ci.up.Menf[i], V1.ci.low.Menf[i]))
  mtext(V1.names.Menf[i], side=1, at=i, 1.12, cex=1.9)
}
points(1.15:3.15, V1.means.Womenf, pch=19, cex=1.5)
for(i in 1:length(V1.means.Womenf)){
  lines(c(i+.15,i+.15), c(V1.ci.up.Womenf[i], V1.ci.low.Womenf[i]))
}
text(2.9, .85, labels = "Men", cex=1.5)
text(3, .75, labels = "Women", cex=1.5)
points(2.7, .85, pch=22, cex=1.5)
points(2.7, .75, pch=19, cex=1.5)
# Export to PDF as 7x4

# Generate in-text quantities fo section 4 of the appendix
Fairness.Fav.v.Unfav <- lm(V1_fair ~ men*V1_fav , data=dataV1_fair[dataV1_fair$V1_equal!=1, ])
summary(Fairness.Fav.v.Unfav) 
# Favorable effect amongst women (0.357, p<0.001)
# Interactione effect of men*Favorable (0.401, 0<0.001)

# For section 5 of the appendix

# Test effect of thinking of fairness as rewarding hard work, as opposed to the baseline of fairness as equality of outcomes OR equality of opportunity (combines both equality measures into the baseline)
dataV1_fair$reward.work <- ifelse(dataV1_fair$fair4==2, 1, 0)

# Table 4 of appendix
mod.fair.work.c <- lm(V1_fair ~ V1_equal*reward.work + V1_fav*reward.work, data=dataV1_fair)
stargazer(mod.fair.work.c) # order of variables is manually adjusted

# Black Respondents:How many selected "rewarding those who work the hardest" 
table(dataV1$fair4[dataV1$race==2 & dataV1$fair4==2]) # 38 

# For Table 5 of appendix
# Compare the distribution of White and Black respondents' beliefs 
# about what fairness means to them, from our full sample, to the White and Black 
# respondents in Brutger and Rathbun's original data


# Data manually entered into appendix Table 5
# White respondents in full sample
table(data$fair4[data$race==1]) # 2277 1316 1283 = 4876

# Black respondents in full sample
table(data$fair4[data$race==2]) # 227  80   127 = 434

2277 /4876 # 47% whites choose "treating everyone equally"
227/434 # 52% Blacks choose "treating everyone equally"

1316/4876 # 27% whites choose "rewarding those who work the hardest"
80/434 # 18% Blacks choose "rewarding those who work the hardest"

1283/4876 # 26% whites choose "Helping those most in need so that they can have the same opportunities as everyone else"
127/434 #29% Blacks choose "Helping those most in need so that they can have the same opportunities as everyone else"


# White respondents in Brutger/Rathbun Sample
table(dataV1$fair4[dataV1$race==1]) # 1176 666 614 = 2456

# Black respondents in Brutger/Rathbun Sample
table(dataV1$fair4[dataV1$race==2]) # 102 38 69 = 209

1176 /2456 # 48% whites choose "treating everyone equally"
102/209 # 49% Blacks choose "treating everyone equally"

666/2456 # 27% whites choose "rewarding those who work the hardest"
38/209 # 18% Blacks choose "rewarding those who work the hardest"

614/2456 # 25% whites choose "Helping those most in need so that they can have the same opportunities as everyone else"
69/209 #33% Blacks choose "Helping those most in need so that they can have the same opportunities as everyone else"


# Before rescaling variable, plot historgram of national attachment for Whites and Blacks
# Appendix Figure 2: export as 4x4 size for each panel
xname <- "National Attachment"
hist(dataV1_fair$NatAtt[dataV1_fair$race==2], breaks = 9, main = "Black Respondents", xlab = xname)
hist(dataV1_fair$NatAtt[dataV1_fair$race==1], breaks = 9, main = "White Respondents", xlab = xname)

# The following lines are for Appendix section 6 discussion
mod.nat.att <- lm(NatAtt ~ white, data = dataV1_fair)
summary(mod.nat.att) # 0.29 p<0.019 (whites have slightly higher NatAtt)
summary(dataV1_fair$NatAtt[dataV1_fair$race==2])
summary(dataV1_fair$NatAtt[dataV1_fair$race==1])


# Rescaling function for more easily interpretable results
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

# Scale cooperative internationalism and  national attachment measures for easy interpretability
dataV1_fair$ci <- rescale(dataV1_fair$ci)
dataV1_fair$NatAtt <- rescale(dataV1_fair$NatAtt)

# Interaction model for Table 6 of appendix
mod.1.n.int.c <- lm(V1_fair ~ V1_equal*NatAtt + V1_fav*NatAtt, data=dataV1_fair)
stargazer(mod.1.n.int.c) # order of variables is manually adjusted in appendix

# Appendix section 7
# Consolidated ideology and partisan models for Table 7 of appendix
dataV1_fair$Democrat <- ifelse(dataV1_fair$party==1, 1, ifelse(dataV1_fair$party==2,0, NA)) #only those who are Democrat or Republicans are coded
Party.Fairness.c <- lm(V1_fair ~ Democrat*V1_equal + Democrat*V1_fav, data=dataV1_fair[is.na(dataV1_fair$Democrat)==FALSE, ])


dataV1_fair$liberal <- ifelse(dataV1_fair$ideology<3, 1, ifelse(dataV1_fair$ideology>3,0, NA)) #only those who are liberal or conservative are coded
Ideology.Fairness.c <- lm(V1_fair ~ liberal*V1_equal + liberal*V1_fav, data=dataV1_fair[is.na(dataV1_fair$liberal)==FALSE, ])

# Table 7 of appendix
stargazer(Party.Fairness.c, Ideology.Fairness.c)

# Section 8 of appendix
# Test whether education could explain differential results
# edu2:1=less than high school or high school / GED; 2=some college, 3=2 or 4 year college degree; 4=masters, doctoral, or professional degree 

# Create indicator for college education
dataV1_fair$college <- ifelse(dataV1_fair$edu2>2, 1, 0)

# Consolidated model for appendix
College.Fairness.c <- lm(V1_fair ~ college*V1_equal + college*V1_fav, data=dataV1_fair)
# Table 8 of appendix
stargazer(College.Fairness.c)


# For section 9 of the appendix
# For marginal effects, first run "interaction_plots.R" to generate functions created by Anton Strezhnev 06/17/2013

# Rename variables for figure
dataV1_fair$Fairness <- dataV1_fair$V1_fair
dataV1_fair$Favorable <- dataV1_fair$V1_fav
dataV1_fair$Unvaforable <- dataV1_fair$V1_unfav
dataV1_fair$White <- dataV1_fair$white

White.Fairness.Fav.v.Unfav <- lm(Fairness ~ White*Favorable, data=dataV1_fair[dataV1_fair$V1_equal!=1, ])

# Appendix Figure 3: export as 4x4
par(mar=c(4,5,2,1))
interaction_plot_binary(White.Fairness.Fav.v.Unfav, effect="Favorable", moderator="White", 
                        interaction="White:Favorable", title="Favorable vs. Unfavorable",
                        xlabel="", ylabel="Marginal effect of Favorable Treatment \n vs. Unfavorable Baseline", factor_labels=c("Black","White"))


# Appendix section 10
# Look at effect of working in import competing industry
# We matched the respondent's employment industry to trade and production data from the Organization of Economic Cooperation and Development (OECD)
# Data is from the OECD 2015 data available at https://stats.oecd.org.
# We then coded an indicator variable, called import-competing, which equals one for sectors with an import-share in the top quartile of the import-shares for all sectors.
dataV1_fair$imp_comp <- ifelse(dataV1_fair$industry>4 & dataV1_fair$industry<9, 1, 0) 

mod.imp1.c <- lm(V1_fair ~ V1_fav*imp_comp + V1_equal*imp_comp, data=dataV1_fair)
# Table 9 of appendix
stargazer(mod.imp1.c) # the order of variables is manually adjusted


# Appendix section 11
# Analysis for Appendix Figure 4

# Create data subset to those who answered the support and fairness dependent variable
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE ) 
# Create data subset those who answered fairness question
dataV1_fair <- subset(dataV1[is.na(dataV1$V1_fair)==FALSE, ])

# Mediation analysis comparing favorable versus unfavorable agreements
# Subset to favorable and unfavorable treatment
set.seed(43216)
dataV1fav.unfav <- subset(dataV1[dataV1$V1_unfav==1 | dataV1$V1_fav==1, ])
dataV1fav.unfav.W <- subset(dataV1fav.unfav[dataV1fav.unfav$white==1, ])
dataV1fav.unfav.B <- subset(dataV1fav.unfav[dataV1fav.unfav$white==0, ])

## Begin with mediation analysis for white respondents
## fairness of agreement as a mediator
fair3.med.W <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.W)
fair3.dv.W <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.W)
fair3.med2.W <- mediate(fair3.med.W, fair3.dv.W, treat="V1_fav", mediator="V1_fair", sims=1500)

## Sensitivity test for left panel of Appendix Figure 4, export as 5x5
sens.fair.med.W <- medsens(fair3.med2.W)
par.orig <- par(mfrow = c(2,2))
plot(sens.fair.med.W, main="", ylim=c(-1,1))

## Now mediation analysis for Black respondents
## fairness of agreement as a mediator
set.seed(43216)
fair3.med.B <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.B)
fair3.dv.B <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.B)
fair3.med2.B <- mediate(fair3.med.B, fair3.dv.B, treat="V1_fav", mediator="V1_fair", sims=1500)

## Sensitivity test for right panel of Appendix Figure 4, export as 5x5
sens.fair.med.B <- medsens(fair3.med2.B)
par.orig <- par(mfrow = c(2,2))
plot(sens.fair.med.B, main="", ylim=c(-1,1))

